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Abstract 

We consider the Euler equation of quasi-geostrophic fluids which is widely used in 
weather forecast. Our goal is to study explicit volume-preserving numerical methods 
for very long simulations on an energy and enstrophy preserving discretization. To 
this purpose, we compute the average fields and estimate statistical parameters. It is 
observed that the statistical parameters depend on the integrals of the discretization 
and that the computed parameters are affected by the error in those, even if the volume 
of the phase space is correctly preserved. We conclude that monitoring the error in 
the integrals by the explicit volume preserving method can be used as an indication 
of how good the estimated parameters are. 

1 Introduction 

This paper deals with the problem of computing long-time trajectories for the Euler equa- 
tions of quasi-geostrophic fluids by means of volume preserving numerical methods. Such 
equations are widely used in applications like weather and climate predictions on Earth 
as well as for other gas planets, such as Jupiter. Due to the highly inhomogeneities in the 
atmosphere or ocean, there is a wide range of unresolved scales of motion and the system 
displays a seemingly chaotic behaviour. As a consequence, it is impossible to simulate 
particular solutions with a prescribed accuracy. The application of ideas from statistical 
mechanics to these geostrophic flows, on the other hand, has lead to promising predictions 
for the ocean, atmosphere and giant planets, in agreement with contemporary observations 
|MW06j . 

The statistics are obtained by performing very long numerical simulations for some 
initial condition. Because of the seemingly chaotic behaviour and the large number of 
unresolved scales, the numerical approximation produced will be far off the trajectory one 
was set to simulate. However, as long as the iteration map is ergodic, time averages are 
equal to space averages almost everywherej^ If the iteration map is ergodic, the long time 
trajectory will visit a large number of points in the phase space with the same probability, 
and the time averages will converge to the spatial averages. Ergodicity is equivalent to the 
preservation of an invariant measure (preservation of a volume form) . Thus it makes sense 

^ Under some continuity assumptions on the map, the measure and the characteristic function, the 
"almost everywhere" convergence can be replaced by "everywhere" convergence. 



to study the effect of volume preserving integrators on the long-time dynamic of such 
systems. Our starting point will be Arakawa's energy and enstrophy preserving semi- 
discretization |Ara66) . which has been shown to preserve volume in phase space [DF07]. 

We will investigate an explicit volume preserving numerical integrator based on split- 
ting methods. The method does not preserve exactly the underlying integrals of the semi- 
discretization, however it is self-adjoint and therefore expected to preserve some nearby 
integrals by virtue of backward error analysis. In order to construct an explicit splitting 
method, we will show that the iih component of the vector field does not depend on the 
ith variable - thus the differential equation is diagonal free in addition to be divergence 
free (the latter was shown in |DF07] ). This will allow us to split the vector field in canon- 
ical shears, which are volume preserving, see |MQ04 MMKQZ08] . To retain the right 
complexity, the splitting must be performed in such a way that computing one step of the 
explicit method does not cost more than evaluating the vector field (0(A^^) operations, 
where N'^ is the number of variables), which is a tricky task due to the presence of the 
Laplacian operator, that mixes up all the variables. This will be discussed in Section [5| 
The simplest way to implement the canonical splitting is by solving in a sequential or- 
der for the variables. However, our preliminary numerical tests showed that this yields 
a poor method in terms of error propagation. It is known that chessboard-like ordering 
of the variables can yield better error propagation for some applications (see for instance 
|FHL97| ). We will also propose a new ordering method (see ^5.2) that is obtained by min- 
imizing a measure of the error of the splitting method: we identify the commuting vector 
fields, collect them together, and add, at each step, those that whose commutation error is 
minimal among the remaining one. This aggregation gives a significant improvement for 
the error accumulation of the numerical simulations. We perform a number of long term 
simulations, for the different indexing strategies and different number of variables in ^ 

The conclusions of this work can be summarized as follows: for a relatively small 
number of discretization points, the second order explicit volume preserving numerical 
method gives satisfactory results when implemented with our indexing strategy, and the 
errors on the conserved quantities stay small over a very long integration time, yielding 
statistically reliable estimated parameters. However, when the number of discretization 
points increase, there is a deterioration of the long-time conservative properties of the 
methods, even with our proposed indexing strategy. Thus, it might be appropriate to 
increase the order of the method. Although it is disappointing to observe this effect of the 
dimensionality on the long-time conservative properties of the proposed volume preserving 
scheme, there is also a constructive side, as the error on the conserved quantities can be 
used to monitor the error on the statistical quantities of interest and hence is a measure 
of their reliability. 



2 The quasi-geostrophic model 

We will consider one of the simplest geophysical models, a barotropic two-dimensional 
fiow with topography on a periodic domain. The model is described by the differential 
equations 

q = A^p + h, |f+V'0xVg = O, 
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where q is the vorticity, V is the stream function and Vu = [dxU, dyv)^ and h is the 
topography. All functions are defined on a a square [0, 2tt] x [0, 27r], and we assume periodic 
boundary conditions. Due to the commuting properties of the differential operators and 
the differentiation rules, the equation can be written in the form 

§-:^A^('?,V') = o, (1) 

where Jm is any of the following equivalent formulations 

Jo{u,v) = VuxVv (2) 
j£{u,v) = Vx{uVv) (3) 
Jziu,v) = -Vx{vVu). (4) 

The system is Hamiltonian, with Poisson bracket 

{F,G} = j qJiVF,VG)dx (5) 
and Hamiltonian energy functional 

The Poisson bracket is degenerate and the system possesses an infinite set of Casimir 
invariants, of the form 

J /(g) dx, 

where / is any scalar function. It is customary to consider invariants of the form J dx, 
for n = 1, 2, . . . (vorticity moments). The most important are the the first moment, also 
referred as the total circulation 

C{q)= [ gdx (6) 



(or total vorticity), and the second moment of vorticity, also called enstrophy, 

Z{q) = \jq'd^. (7) 

The higher moments do not seem to have much relevance in applications, except for the 
third moment 

o / g^dx. 
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which has been seen to have some relevance in the statistical analysis of the discretized 
flow |AMn3] . 
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3 The conservative Arakawa decompositions and volume 
preservation 

Arakawa |Ara66| considered semi-discretizations of the above equations, in which the dif- 
ferential operators were replaced by standard central divided difference operators: 

dx ~ Dx, dy ^ Dy. 

In the discussion that follows, it is important that the operators are skew-symmetric, 
i.e. = —Dx,Dy = —Dy. We consider a symmetric discrete approximation to the 
Laplacian, A^r = A^, and, for the time being, we assume it to be invertiblej^ We introduce 
a discretization of the domain. Thus, for an arbitrary function u, we set Uij ~ u{i5,j5), 
where 5 = 2tt/N is the spacing in the x and y directions (isotropic discretization) and 
N is the number of mesh points in each direction (for a total of A^^ grid points and as 
many variables). Further, by u, we will denote the vector of the Uij using the standard 
column-wise ordering. Letting 

V' = A^l(q-h), 

we can express the equations in terms of the sole variable q, so that the semi-discretized 
system takes the form 

^=Jm(.ci), (8) 
where Jm is one of the expressions below: 

Jo(q) = (D.q)*(I),A^l(q-h))-(I),q)*(Z),,A^l(q-h)), (9) 
Jsiq) = L>,(q*Dj,A^l(q-h))-L>^(q*Z),A^l(q-h)), (10) 
Jziq) = Dy{{DxCi)*A-^\q-h))-D.,{{DycD*A-^\ci-h)), (11) 

A4 = 0,£,Z, which are finite difference discretizations of the corresponding formulations 
Q-Q. The '*' operator denotes componentwise multiplication of vectors, (u*v)j = UiVi, 
so that the product 

u^(v * w) = UiViWi = v"^(u * w) = w'^(u * v) 



I 



is fully symmetric. Differently from the analytical case, the semi-discretizations ([9])-(ll) 
are not equivalent, rather, they have very different properties. All of them preserve the 
discrete version of the total vorticity, 

C(q) = q'^152, (12) 

otherwise, the flows have different conservation properties. Arakawa showed that the 
choice ( 10 ) preserves a discrete version of the energy, 

i?(q) = -^^^(q-h)<52, (13) 



^When this is not the case, a symmetric pseudo-inverse A^^ wiU be used instead, see discussion later. 
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(but not the enstrophy), while (11) preserves the discrete enstrophy, 



^(q) = 



(14) 



(but not the energy). Jq preserves neither the discrete energy (13) nor the enstrophy (14). 



In addition to ([9])-(ll), Arakawa considered also a semi-discretized flow that preserves 
energy and enstrophy (and also linear momentum), 



Jezici) = ^(Jo(q) + (q) + JzicD). 



(15) 



The proofs of conservation can be found in the original work of Arakawa. A more modern 
approach, based on Nambu bracket formalism, can be found in |DP07j. 



The semi-discretizations (^9^-(ll) and (15) are also volume preserving. A proof, using 



the skew-symmetry of the discretized differential operators , Dy and the symmetry of 
Ajf , as well as the properties of the trace operators, can be found in |DF07] . 



4 Volume preserving numerical methods 

The divergence-free differential equation 

q = f(q), v^f = V^=0 

q(0) = qo, ^dqi ^' 

preserves the volume, as it has the property that det9q/9qo = 1. Similarly, a numerical 
method for q = f(q), 

q„+i = *^,f(qn), n = 0,l,..., 
is volume preserving provided that 

det^^ = l. (16) 

or, equivalently, det[9q„+i/9q„] = 1 for all n (whenever V-^f = 0). Following |MMKQZ08| , 
a system obeying the condition 

dfi 

-^ = 0, z = l,2,..., 

is called diagonal-free system, since the i-th component does not depend on the i-th 
variable. Diagonal-free systems are obviously divergence free, but the converse is not true, 
as there exists divergence- free system that are not diagonal free. It is easier to design 
volume preserving numerical methods for diagonal- free systems, while the general case 
with nonzero diagonal contributions is typically more complicated |MMKQZ08| . 

Let us consider a splitting method, based on the splitting of the function f in M "sub 
vector fields" 

f = pi • • • F^'^ 
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each of them also divergence free. Assume that each vector field F'^ is integrated exactly 
by for k = 1, . . . , M. Then, the composition 



*^,f = *1 o . . . o *f 

is also volume preserving and yields an order one method. We will consider the symmetric 
composition, 

*r,f = ^1/2 ° • • • o ' o *f ° ' ° • • • ° *'/2. (17) 

which yields a second order method (each is an exact integrator, thus self-adjoint, 

A simple and natural splitting of diagonal-free systems in vector fields that can be 
solved exactly is that of canonical shears, advancing one coordinate at a time, while 
keeping the others constant. In detail, we take M to be the number of variables, and set 
Fj to be the vector field equal to for the i-th variable, and for the remaining: for 
i = 1,...,M, 

qj=0, i = l,...,M, j^i. 

It is immediate to see that the variables qj, j ^ i are constant, while the i-th variable can 
be integrated exactly by a step of the Forward Euler method, as fi depends solely on the 
other variables. In other words, the exact solution is given by 

z = *;(q)=e-N: ^^ = 'l^ + rfM), 

Zj =qj, J = 1, . . . ,M, J / I. 

Since solves q = Fj(q) exactly and Fj is divergence free, it is volume preserving. 



Alternatively, one can proof directly the condition (16) by observing that the matrix 
'^^q^^ differs from the identity matrix only for the off-diagonal elements in row i, and 



clearly it has determinant one. Consequently, the symmetric composition (17) yields an 
explicit, symmetric, order two, volume preserving integrator. 

It must be noted that composition methods do not necessarily preserve first integrals 
and/or conserved quantities. A first integral I{t) will be preserved if and only if each 



of the flows also preserves I{t), which is typically not the case of the shears. So, 



the symmetric composition method (17) is not expected to preserve the integrals and the 



Casimirs of the £Z discretization. However, due to the symmetry of the composition (17), 
there exist theorems that guarantee that the conserved quantities are approximately well 
preserved over long simulation times |HLW02l ILR04] . 



5 Construction of the VP numerical method 

In the context of the QG equations, we have f(q) = Jm{^}j so that fi is the ith component 
of Jm- Since 1/7 = A^^(q — h), the Jm vector fields are all quadratic in the q variable. 
This means that they can be written in the form 

Mci) = ^ali{qk-hk)qi = {ci-hfA'q, i=l,...,N\ (18) 
k,l 
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where = (a^^)- We have omitted the dependence of the matrices and their elements 
on the vector field Ai = 0,£,Z, £Z for sake of simiplicity. In order to be able to apply 
the canonical shears methods described above, we must make sure that the vector field is 
diagonal-free. To do so, we study the coefficients a^^. 

The coefficients a^. ^ can be computed explicitly and once for all at the beginning of 
the numerical simulations. They do not depend on the topography or the vorticity, but 



only on the differential operators. Let Lx = DxA^^ and L. 



Then: 



SZ : 



all 



-'xji.kl 



■^{iDx)i,i[{Ly)i^k + iLy)i,k] - iDy)i^i[{Lx)i,k + {Lx 

+ [{Dx\:*{Dy)l,-{Dy\,*{Dx)l,-\ (A^^):,fc ), 



row vector 



col. vector 



: 


ah 


= (Dx) 


i,l{Ly)i,k - {Dy)i i{Lx) 


i,k 


S : 


a\,i 


= (Dx) 


i,l{Ly)l^k — {Dy)i,l{Lx) 


l,k 


Z : 


4,1 


= [{Dx 


)i, * {Dy)l, - (Dy),, * 


{Dx 








row vector 





I (A^)^ > 
col. vector 



(19) 

(20) 
(21) 
(22) 



where we have used Matlab semicolon notation to denote rows/columns of matrices. 

We use the standard five point stencil with periodic boundary conditions. Let E the 
N X N shift matrix (periodic). Then, 



D 



E-E^ 
26 ' 



where 5 = 2'k/N. Similarly, 

E-21 + E'^ 



D2 = 



5^ 



Ajv = -D2<X)/ + /(8)-D2• 



(23) 



For periodic boundary conditions, the matrix D2 is circulant hence singular. In turn, 
A AT is block circulant, with circulant block, and singular. Further, if the matrix D2 is 
symmetric, also the blocks in Ajv are symmetric. 

The matrix A^v has a rank-one subspace corresponding to constant solutions. Thus, we 
need an "inverse" that projects constant functions to zero. The Moore-Penrose pseudo- 
inverse matrix A^ has this property. Let 



A 



N 



VAV^ 



the an orthogonal eigenvalue decomposition of A^v (which exists, since Ajv is symmetric, 
hence normal). One eigenvalue of A^v, hence one diagonal element of A, will be equal to 
zero. Define the matrix A+ with diagonal elements 



1/Ai 




The pseudoinverse is then defined as A^ 
Some important observations are in place 



i if K,i 7^ 
otherwise. 

VA+V^. 
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Proposition 1 Consider the Jm vector field in (l8|), whose i-th component coefficients are 



given as in (18). The matrices, obtained using the divided differences (23) operators, are 
sparse and have only a few number of columns, depending on the choice of the discretized 
differential operators. In particular, for the 0, £, Z case, only 4 nonzero columns, while 
for the EZ case, 8 nonzero columns. 

Proof. This can be proved as follows. For a fixed i, it is readily observed that {Dx)i^i is 
nonzero for two values of / (and independently of A;), and so is {Dy)i i, yielding possibly 
nonzero elements for four columns for the 0,£" cases. Further, the row vector \{Dx)i,: * 
{Dy)i . — {Dy)i^- * (D-c); .] can be nonzero at at most four values of I (two values from the 
first product, further two from the second product). This implies that the Z-case matrices 
also have only four columns. Finally, for the £Z cases, we have the superposition of the 
0,1? case and the Z case, giving 8 columns (at most). I 

Proposition 2 /i(q) does not depend on qi, i.e. a* ^ = ^ = for all k,l = 1, . . . , N'^ . 

Proof. In the case : for / = i, it is obvious that • = since P)x i T^y are skew 

symmetric hence they have zero diagonal elements. In the case k = i, we have {Lx)i^i = 

{Ly)i^i = 0. We show the case {Lx)i^i = 0, the other one is similar. We have (L^ 
^Tt — ^Tn A-i^. — ^^Tn A-l'^.^T _ ^TA-in ^. _ ^Tn a-i 



The third last passage follows from the symmetry of Atv and skew-symmetry of Dx and 
the second last follows from the fact that the differential operators Dx and Ajy commute. 
Hence {Lx)i,i = -{Lx)i,i = 0. 

In the £ case : for / = i, this is the analogous of the case. For k = i, we have ^ = 

{Dx)i,l{Ly)l^i - {Dy)i^i{Lx)l^i = {Dx)i,l Y.mi^y)l,m{^'^)m,i - {Dy)i^l ^rni^^)l,mi^N^)m,i- 

Because of the tensor-product structure of Dx and Dy, it is immediate to observe that 
{Dx)i,i{Dy)i^l = for any i, I. Hence at most one of them is nonzero at each time. 

Assume {Dx)i,i / (hence {Dy)i^i = 0). We look at ^rn(^y)i,m{^J^)m,i- Because of 
the structure of Dy, for any fixed index I, there are at most two indices mi,m2 for which 
{Dy)i^rni and {Dy)i^m2 are nonzero. Moreover, {Dy)i^mi = —{Dy)i^m2-, hence it is sufficient 
to show that [/\'^)mx,i = (A^^)m,2,i- This latter statement follows because of the structure 
of A]^^, which has circulant and symmetric blocks (in addition to being block-circulant 
and block symmetric), since such is A^r. The case {Dy)i^i ^ (hence {Dx)i,i = 0) can be 
proven in a very similar manner. 



D^A-^ei = {ejDx^-^^eiY = -ejAj^Dxei = -ejDxA-^Bi = -{L^ 



In the Z case : for / = i, the row vector in (22) is obviously identically zero, since, as 
observed above, {Dx)i,m{Dy)i^rn = for any i,m. 

For k = i, the proof is a verification as in the £ case. For any i, I, either the row vector 
is identically zero (for instance if the indices i, I correspond to the same block), or, fixed i, 
there exist (at most four) values of / for which the row vector {Dx)i^: * {Dy)i^- has a single 
nonzero element, corresponding to the index mi. For the same values of i,l, the term 
(Dy)i^- * {Dx)i^: has exactly the same single nonzero element, but at a position m2. It is 
true that (A^v)"^ ■ = (A^)^^ j. The proof is tedious and technical, but, again, it follows 
from the special structure of (Ajv)~^, being block circulant, with circulant blocks, block 
symmetric, with symmetric blocks, and its connections with convolutions. 

Finally, the £Z case is the combination of the three cases, hence it follows from those. 
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Proposition 3 The matrices in (18) are the identical up to a shift on the physical 
domain. 



Proof. The coefficient of the matrices are determined only by the relative position of qk, 
qi with respect to the variable qi. This is true because the differential operators in the 
definition of the vector field: the interaction between q^ and qt depends only on their 
position in the grid. I 

5.1 Complexity of the VP scheme for the 8Z discretization 

From this point, we focus on the £Z case (energy and entropy-preserving discretization, 



whose j4* matrix elements are given in ( 19 ). The matrices (in fact, because of Prop (|3|), 
they can all be identified with a single matrix, that is then shifted over the domain) 
as well as the differential operators are computed once and for all. Each evaluation 
requires a vector-matrix-vector multipication, where the matrices have only 8 nonzero 
columns and N"^ rows. Hence each fi requires 8A^'^ multiplications and as many additions. 
This must be repeated for all the A^^ variables, and the result doubled if a symmetric 
composition is used. Hence the complexity of one step of the symmetric method is 32 A^^, 
counting both multiplications and additions. Thus, evaluating one step of the explicit 
splitting method has the same order of magnitude of computing the whole vector field f . 

5.2 Indexing strategies for the £Z discretization 

Given a splitting of a vector field in M terms, there are Ml possible orderings (permuta- 



tions of indices) for the symmetric composition (17). A standard approach would be that 
of choosing a sequential ordering. Other popular approaches are based on checkerboard 
orderings. Our preliminary numerical simulations indicated that these strategies were not 
satisfactory, hence a better strategy for the ordering of the vector fields was needed, as 
each ordering will have a different effect insofar error propagation is concerned. 



To have an idea of the error propagation by the splitting method, we analyze (17) by 
the symmetric BCH formula. 



e^/2Fl . . . gT/2FMgr/2F„ . . . gr/2Fi ^ ^Z(t) 



where 



M 3 M Z ^ 

Z(r) = rj;F,-^ [F,, [F,,Fd] - ^ [F„[F„F,]] + O(r^). 

*=1 i,j,k = 1 i,k = 1 

j<i,j<k i<k 

The leading error term is governed by the double commutators of vector fields, defined as 
[F, G] = VF G — VG F, where VF is the Jacobian of the vector field F. It is immediate 
to observe that collecting together commuting vector fields will decrease the number of 
error terms, it is reasonable to collect together commuting vector fields. When the vector 
fields do not commute, minimizing the commutator is still a reasonable strategy. Thus we 
estimate the norm of the commutators and collect together terms in order of increasing 
non-commutativity. To this goal, after computing the matrices A^, for i = 1, . . . , A^^, for 
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each A^, we look at the column values and the row values. If, for an index k = 1, . . . , N'^, 
A'^{:, k) = and A^{k, :) = 0, then /j does not depend on g^,. Similarly, if A^{:, i) = and 
A''{i, :) = 0, then fk does not depend on qi, hence Fj = fief and = commute. 

Disregarding the dependence on the topography h, the vector fields are of the form 
Fj = (q"^^*q)ej (i.e. they are quadratic and advance one variable at a time, qi), hence 
Fj and Fj will nearly commute if the coefficients in A^ that invoke the j-th variable are 
small. Thus, we look at 

fc=i 

The variable c*- is a measure of the independence of the vector field Fj on the variable qj. 
More formally, we have 

■T 

VF, F, = (q^A^q)e,q^^^^±^e„ 

hence, assuming ||q||2 < 1 and 11^4-' II2 < a, we can bound 

ll[F^,F,]||2<a4 

for any i, j = 1, . . . , N. If Aj = A\j = for all k, then Fj is independent of qj and Fj 
and Fj commute. 

The vector c* = [c^, . . . ,c]y2]"^ is then reshaped to the grid, giving a matrix C*. The 
matrix C* is essentially the same for all the variables, as the matrices A^ are essentially the 
same for all the variables (cfr. Prop.js]), they only differ by a shift in the location. Boundary 
wrapping must be considered, as we the problem has periodic boundary conditions, see 
Figure [3} 

For simplicity, we illustrate the case when N = 2^. The arguments extend also to the 
generic case even (our numerical tests are all performed with even) . The C* matrices 
look like perfectly symmetrical 'vulcanos'. The locations at the same distance from the 
centre have the same commutation weight. The commutation weight decreases with the 
distance to the centre and it is minimal for the point on the grid that is furthest away 



(considering wrap-around) from it. For the standard divided difference choice (23) we get 
the following commutation pattern: divide the grid in four quadrants, NW, SW, NE, SE. 
The flow corresponding to each variable i in the NW quadrant commutes with the flow 
of the corresponding variable in the SE quadrant. Similarly, the flow of a variable in the 
SW quadrant commutes with the flow of the corresponding variable in the NE quadrant 
(see figures [l]j2] for an illustration in the case N = 8, i.e. 8x8 grid, N'^ = 64 variables). 
Commutation can be proved, but the proof (technical, in the style of Prop. [2]) is not so 
relevant for the discussion. Notice also that the closest positions (hence variables) to the 
centre are always those with the highest non-commutativity. This explains why a plain 
sequential indexing of the variables, 1,2,3,..., N"^, will typically display much larger error. 

Once the matrix is computed, each C* is obtained by shifting on the relevant 
position in the domain. If we consider the indexing {l,i), the error will be approximately 
the sum of the two commutation weights, i.e. Cust = + (7*. Since we aim at minimizing 
these cumulative sums, we choose the variable 12 to be the position where has a 
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Figure 1: The commutation diagram NW-SE (left) and SW-NE (right) for the vector fields 
Fj = fi^i for the choice of the standard 5-points discretization of the Laplacian by divided 
differences. 



minimum value (different from the centre of the vulcano). This gives us a partial ordering 
{h, ^2) = (1) ^2); with cumulative commutation weight Cust = + C*^. Then we look at 
the values of Ciist, sorted in increasing order, ignoring the minima corresponding to 1,^2; 
and take the next minimal value. In case of a tie, we obtain a new list, of candidate indices 
for which we compute the cumulative weight (starting from the first index in the new list). 
The new list, sorted according to the cumulative weight, is then added to the previous 
list, and the corresponding commutation weight matrices are added to Cust- We look at 
the new minima, get a new list of candidate indices, eventually sort them, and so on, until 
the final ordering list is obtained. The procedure can be formalized as an algorithm. 



Algorithm 1 [MinCom] Indexing of the vector fields for divided differences (23), by min- 
imal cumulative estimated commutation weight. 

Description: Sorting of the vector fields. Creates a ordered list of indices 
from 1 to N"^ (number of vector fields). 

Input: The commutation weight matrix , N, starting index il. 

Output: A list of indices (1,^2) • • • j^atz) with a minimal cumulative estimated commutation weight. 

% Initialization 
list =[il]; 
newlist = list; 

Clist = zeros (N ,N) ; % "cumulative commutation" weight matrix 
% Main loop 

until all variable {1, . . . , A^^} are in list do 
for all i in newlist 

Shift the commutation weight matrices in position i to obtain C*; 
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Figure 2: A matrix showing the commuting vector fields corresponding the standard 5- 
points discretization of the Laplacian by divided differences. The entry in the position 
is 1 if fi does not depend on qj, zero otherwise. This graph displays the same 
information as the NW-SE and SW-NE diagrams in Figure [TJ 

Clist = Qist + C**;' 

end for 

% Find the next most commuting indices newlist (candidate indices) by getting 
% new minimal value: 

[sval, newlist] = sortCCustC: )); % sval temporary variable with sorted values 
Remove from newlist all indices already in list; 
Set inewlist to be the first index in newlist; 
minval = sval (inewlist) ; % new minimum 

Find all indices in newlist that have value equal to minval and remove rest; 
Sort the element in newlist so that they locally . . . 

have the least cumulative commutation matrix with starting . . . 

index inewlist (recursive call) ; 
list = [list, newlist] ; % Concatenate lists 
end do 

The procedure for the algorithm and the intermediate cumulative weight matrices for the 
N = 8 case is illustrated in Figure |4] below. Interestingly enought, the results of this algo- 
rithm returns a grid that reminds of multigrid approaches, which might give an intuition 
of why this ordering should give better results. But there are differences: in the multigrid 
approach, the computed values are not variable values, rather, they are transferred to the 
relevant variables at the higher resolution by some interpolation operator, until the finest 
level. In this approach, we do compute the actual variable values and no interpolation is 
taking place: all the computations are done at the same grid level. 
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Figure 3: Illustration of the matrices C* with divided differences as in (23), for i = 5,37, 
for a 8 X 8 grid. The values are the same, just shifted across the domain, centered on the 
grid position of the ith variable. F5 commutes with variable F33 and F37 commutes with 
Fi (cfr. Fig.[l]). Notice the wrapping along the boundaries. 



6 Numerical results 

We perform numerical simulations with the J^z vector field on the test problem of |AM031 
IDF07| . The topography function is 

= 0.2COSX + 0.4cos2x. (24) 

The coefficients of our VP implementation based on splitting methods are given in ( |19[ ). 

We test different sizes of the grid. For each grid size, we use the same initial condition, 
computed numerically using an optimization technique and random initial condition, with 
the constraints E = 7, Z = 20, and ffist and third moment equal to zero as in |DF07j . 

We estimate the mean fields of the numerical simulations and compare with the mean 
field predictions of the £Z semi-discretization. 

The numerical experiments consist in estimating the average fiels (1/') and the corre- 
lation /i between the average vorticity (q) and (t/?). In facts, under the assumptions of 
energy and enstrophy conservation, and conservation of linear momentum (total vorticity), 
the statistical theory asserts that 

(q> = /^(^>, (25) 

and, further, that is a function of x only (as such is h{x, y)). The least squares solution 
for gives the estimate 

For the values of energy and enstrophy above, E = 7, Z = 20, and zero total vorticity, 
the predicted value of // was estimated by |DF07j to be ^ « —0.730 for = 22. We have 
computed the predictions for different values of N below. The algorithm for the prediction, 
described in |DF07j and adapted from the spectral analysis in [CF87) , is reported here for 
convenience. The ensemble averages for energy and estrophy are split into a mean part 
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Figure 4: An illustration of Algorithm [T] for indexing of a 8 x 8 grid, showing the cumulative 
commutation matrix Cust with starting index ii = 1. Blue values: commuting/almost 
commuting. Red values: high dependence on these variables (the commutator is expected 
to be larger). Top left to bottom right: At the fist step, newlist=[37] is added (commuting 
with 1). At the second step, newlist = [5, 33]. As they too commute in between, no sorting 
is needed. At the third step, newlist = [3, 39, 7, 35, 21, 49, 53, 17]. They have been sorted 
so that the cumulative commutation matrix of the new list, starting with the inewiist = 3, 
has the smallest possible value in the corresponding positions. At step 4, newlist =[19, 
55, 51, 23], (also sorted). At step 5, newlist = [10, 46, 42, 14, 44, 16, 48, 12, 28, 64, 60, 
32, 30, 58, 26, 62] (sorted) and finally at step 6, newlist =[2, 38, 34, 6, 8, 36, 40, 4, 56, 20, 
52, 24, 22, 50, 18, 54, 29, 57, 61, 25, 9, 45, 13, 41, 47, 11, 43, 15, 27, 63, 31, 59] (sorted). 



and a fluctuation, (En) = Ej\f{{q)) + E'j^, {Zn) = ^^{{q)) + Z'^. Then using a spectral 
truncation, one has 

fc,Z=-Ar/2+l ' k,l=-N/2+l ^ 

N/2 r,,- ,9 N/2 r, 

= 2 2. (mTFT/¥ ' ^ = 2a 2. ^ + P + 

k,l=-N/2+l ' k,l=-N/2+l ^ 

Starting with suitable guesses for /i, a, the actual values are computed by iteration (non- 
linear least squares) using (p6|)-(27) in conjunction with the assumption that E"^ = E = 7 
and Z*^ = Z = 20, see Table|l} 

To compute our averages, we perform a long time numerical simulation. Assuming 
ergodicity (because of volume preservation of the numerical method), we compute the 
expectations relative to the microcanonical ensamble. On the other hand, the theory 
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N 


6 


8 


10 


16 


22 


32 


64 


predicted fi 


-0.3995 


-0.5646 


-0.6402 


-0.7106 


-0.7298 


-0.7409 


-0.7487 



Table 1: Predicted values of ^ as a function of N for fixed E = 7, Z = 20 and topography 



h{x, y) as in (24) 



predictions determine the expectations for the canonical ensamble. Only in the limit, 
N — )• oo, the averages (hence the expectations) coincide. 

The start for averaging time is Tq = 10^, to allow decorrelation of the initial conditions. 
The averages are then computed at T = 10^, 10^, 10^, according to the formula 

1 r 

- -to JTo 

Numerically, denote by i/?^ the numerical approximation at time t = kr. Assume T = 
Kt, Tq = Kqt. Then, 



1 ^ 



k=Ko+l 

A similar formula is used for the q variable. 

In the experiments, we will consider the following types of (forward) ordering of the 
variables 13 

Plain Denotes the standard sequential ordering 1,2,3,... 

BW Checkerboard approach. For some types of problems, this ordering might give sepa- 
ration of variables and better error propagation. In this problem, the BW approach 
does not give variable separation, due to the "mixing" effect of A^. 

MinCom The indexing obtained by the Algorithm [T] computing the minimal cumulative 
estimated commutation weight. 

6.1 Experiments with 8x8 grid 

In the case of the 'Plain' symmetric splitting, we observe a linear decrease of the relative 
energy over long integration time. At T = 10^, the relative energy error is of the order of 
10~^, and this is reflected in the estimate of the //, that also increases at the same pace. 
The relative enstrophy error grows at a smaller rate, but still has a linear behavior. The 
absolute total vorticity (first moment) error also grows, at a much faster pace. 

As noted in |DF07j . a nonzero total vorticity ^ ■ ^ qij will give a constant displacement 



in the relation (25). The total vorticity is a linear integral, and is automatically preserved 
by any linear integration methods. However, this is not the case of splitting methods, 
unless the integral is preserved by any of the split fields. 



^The actual indexing is a symmetric version of the above, i.e. forward and backward, in order to obtain 
an order 2, symmetric, splitting method. 
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Figure 5: 8 x 8 grid, with 'Plain' ordering of the indices, 1,2, ... , N"^ , N"^ — 1, ... ,2,1. Top 
left: Estimated /i. Top right: Average field, and linear fitting plots, {q) vs. {tjj). Bottom 
left: relative error in the energy (solid line) and enstrophy (dashed line). Bottom right: 
absolute errors in the total vorticity (solid line), which is also the first moment of the 
system, and third moment (dashed line). 



We also plot the absolute error in the third moment. This is not as well preserved, 
but bounded. The result is not unreasonable as the semi-discrete equations (15) do not 
preserve the third moment, while first moment, energy and enstrophy are preserved. 

The checkerboard 'BW splitting is shown in Figure [6j We observe that the energy 
and enstrophy are already better preserved than the 'Plain' case. The relative error in the 
energy and enstrophy is still of the order of 10~^ even after lO"^ integration steps. However, 
the total vorticity still shows a noticeable error and decreases to ~ —2 at the end of the 
integration. A comparison of the first moment for the 'BW' splitting and the 'MinCom' 
approach is shown in Figure [8j The absolute error in the third moment is comparable to 
the case of the 'Plain' discretization and appears to be bounded. 

In the case of the 'MinCom' symmetric splitting, we observe a much better and slower 
error propagation over the long time integration. The relative error in the energy and 
enstrophy is still of the order of 10"'^ even after 10^ integration steps. This reflects on the 
estimate of the /U, that does not show a liner increase, as in the 'Plain' case. The absolute 
total vorticity (first moment) error also is much better conserved. The absolute error in 
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Figure 6: 8 x 8 grid, with 'BW ordering of the indices. Top left: Estimated /i. Top 
right: Average field, and linear fitting plots, {q) vs. Bottom left: relative error in 

the energy (solid line) and enstrophy (dashed line). Bottom right: absolute errors in the 
total vorticity (solid line), which is also the first moment of the system, and third moment 
(dashed line). 
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Figure 7: 8 x 8 grid, with 'MinCom' ordering of the indices, as described in Algorithm [T} 
Top left: Estimated //. Top right: Average field, and linear fitting plots, {q) vs. ("0)- 
Bottom left: relative error in the energy (solid line) and enstrophy (dashed line). Bottom 
right: absolute errors in the total vorticity (solid line), which is also the first moment of 
the system, and third moment (dashed line). 

the third moment is comparable to the case of the 'Plain' discretization and appears to 
be bounded. 

6.2 Experiments with 16 x 16 grid 

We repeated the same experiment on a finer grid, N = 16. Again, we tested a 'Plain' 
indexing, a 'BW and a 'MinCom' indexing. The initial solution was the same for the 
three trials. Once again, the BW method failed to provide any reasonable result. The 
splitting using the 'Plain' indexing provided results up to t = 8 x 10^ (blow-up time), but 
the solution exploded thereafter. The estimated fj, shows a clear linear trend until the 
blow-up time. The linear trend is likely a consequence of a linear increase in energy and 
enstrophy by the numerical simulation. 

Finally, Figure [Tl] shows the results of the simulation using the MinCom indexing 
approach in Algorithm [T] The relative errors in the energy and enstrophy, as well as linear 
momentum, are very well preserved, and this reflects in the estimation of See Table [3] 
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Figure 8:8x8 grid. A comparison on the conserved quantities, energy, enstrophy and total 



vorticity of the £Z semi-discretized flow (15) for the 'BW (left) and 'MinCom' (right) 
indexing. 



8x8 


'Plain' 


'BW 


'MinCom' 


T= 10^ 


-0.4283 


-0.4334 


-0.4344 


T = 10^ 


-0.4225 


-0.4294 


-0.4267 


T = 10^ 


-0.3919 


-0.4275 


-0.4262 



Table 2: Estimated values of /i for the different indexing, 8x8 grid. Initial averaging 
time: Tq = 10^. Stepsize of integration r = 0.1. 

for the actual figures. 

6.3 Experiment with 22 x 22 grid 

The same experiment is repeated on a finer grid, N = 22. We tested only the 'MinCom' 
indexing, given the deterioration of the quality of the solution for the two other orderings. 
The initial solution was randomly generated. For this number of mesh points, the energy 
and the enstrophy grow also for this indexing of the variables. The estimated ^ shows a 
clear linear trend. The linear trend is likely a consequence of a linear increase in energy 
and enstrophy by the numerical simulation. See Table |4] for the actual figures. 



16 X 16 


'Plain' 


'BW 


'MinCom' 


T = 10^ 


-0.7023 


-0.6951 


-0.6939 


T = 10^ 


-0.7261 


-0.7061 


-0.6947 


T = 10^ 


NaN 


-0.8013 


-0.6954 



Table 3: Estimated values of /x for the different indexing, 16 x 16 grid. Initial averaging 
time: Tq = 10^. Stepsize of integration r = 0.1. 



19 




Figure 9: 16 x 16 grid, with 'Plain' ordering of the indices. Top left: Estimated /i. Top 
right: Average field, and linear fitting plots, {q) vs. Bottom left: relative error in 

the energy (solid line) and enstrophy (dashed line). Bottom right: absolute errors in the 
total vorticity (solid line), which is also the first moment of the system, and third moment 
(dashed line). 



22 X 22 


'Plain' 


'BW 


'MinCom' 


'MinCom4' 


T = 10^ 






-0.7248 


-0.7245 


T = 10^ 






-0.7342 


-0.7234 


T = 10^ 






-0.8194 


-0.7201 



Table 4: Estimated values of /x for the different indexing, 22 x 22 grid. Initial averaging 
time: Tq = 10'^. Stepsize of integration r = 0.1. The last column is obtained with an 
order 4 implementation of the method, using Yoshida's technique in |6.4[ 
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Figure 10: 16 x 16 grid, with 'BW ordering of the indices. Top left: Estimated /i. Top 
right: Average field, and linear fitting plots, {q) vs. (tp). Bottom left: relative error in 
the energy (solid line) and enstrophy (dashed line). Bottom right: absolute errors in the 
total vorticity (solid line), which is also the first moment of the system, and third moment 
(dashed line). 



21 




Figure 11: 16 x 16 grid, with 'MinCom' ordering of the indices, as described in Algorithm[T| 
Top left: Estimated /x. Top right: Average field, and linear fitting plots, (g) vs. 
Bottom left: relative error in the energy (solid line) and enstrophy (dashed line). Bottom 
right: absolute errors in the total vorticity (solid line), which is also the first moment of 
the system, and third moment (dashed line). 
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Figure 12: Energy, enstrophy and total vorticity for the 'BW (left) and 'MinCom' (right) 
indexing for the 16 x 16 grid. 
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Figure 13: 22 x 22 grid, with 'MinCom' ordering of the indices, as described in AIgorithm[T| 
Top left: Estimated /x. Top right: Average field, and linear fitting plots, {q) vs. (-0)- 
Bottom left: relative error in the energy (solid line) and enstrophy (dashed line). Bottom 
right: absolute errors in the total vorticity (solid line), which is also the first moment of 
the system, and third moment (dashed line). 
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6.4 Simulations using higher order methods 



Clearly, the 22 x 22 simulation, even with the optimized ordering of the splitting terms, 
does not show satisfactory results, as there is a clear linear drift in the estimation of /j,. 
This seems to be apparently due to worse accumulation of error, that eventually displays 
an exponential growth. This is disappointing, as dimensionality seems to play a role in 
how long the method is close to a "nearby problem" . If the energy and the enstrophy (and 
linear momentum) grow too much, it is clear that the computed statistical averages will 
be affected. 

To enforce better preservation of the conserved quantities, we test the approach using 
an order four implementation of the volume preserving splitting method using Yoshida's 
technique |Yos90] : given a self-adjoint method S'^ of order two, it is possible to increase 
the order from two to four using the composition 

= S'q^ o s"^^ o s'^^ 

where a, j3 obey the conditions 

2a + /3 = 1, 2a^ + /S^ = 0. 

There is a unique real solution to the above equations, given by a = 1/(2 - 2^1'^) and 
j3 = —2^/^/{2 — 2-*^/^). By a similar token, it is possible to increase the order from four to 
eight and further, though it is known that this approach might lead to higher leading error 
term and compositions with fewer terms and smaller leading error can be constructed, see 
for instance |BCM08j . 

The results of the order four implementation of the volume preserving splitting method 
for the 22 x 22 example with the same initial condition as the order two method, are shown 
in Figure 14 , and the computed values for /x are tabulated in Table |4j The fourth order 



method has a better conservation of all the integrals of the £Z semi-discretization, and in 
addition it displays a good preservation of the third order momentum. However, there is 
still a mild linear tendency in the estimated ^. 



6.5 Some comments and concluding remarks 

With the increase of the grid dimension, we have observed that the splitting methods 
display an increase in error accumulation, in particular for the total vorticity (first mo- 
mentum), which is an integral of the semi-discrete problem, as are the energy and the 
enstrophy (in contrast to the third moment, which was not an integral of the discretiza- 
tion) . The total vorticity had overall a consistent error propagation worse than the energy 
and the enstrophy. 

The statistical predictions for the QG flow are valid under the assumption that energy, 
enstrophy and total linear momentum are conserved. When this is not any longer the 
case, the numerical computed asymptotic values do not reflect any longer the theoretical 
predictions of the model, and the computed ^ parameter shows a clear linear trend instead 
of converging to an asymptotic value. Thus, particular care has to be used when using a 
splitting method. The ordering of the vector flelds matters for error propagation and we 
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Figure 14: 22 x 22 grid, order 4 implementation with 'MinCom' ordering of the indices, 
using Yoshida's technique. Bottom left: relative error in the energy (solid line) and en- 
strophy (dashed line). Bottom right: absolute errors in the total vorticity (solid line), 
which is also the first moment of the system, and third moment (dashed line). 
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Figure 15: 22 x 22 grid, order 4 implementation with 'MinCom' ordering of the indices. 
Energy, enstrophy and total vorticity. 
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have proposed an algorithm that produces an ordering of the vector fields that is nearly 
optimal. 

When a volume preserving method is used, with nearby-conservation of first integrals, 
the ergodicity property of the flow is respected and the departure from the nearby integrals 
can be used as a criterium of how good the averaged result is, in these case for which the 
statistical quantities depend on the first integrals of the system. If the volume is not 
preserved, the preservation of integrals alone might give false positives, as generally no 
estimation of how much the (high-dimensional) space shrinks or expands is available, 
hence no indication of how reliable the estimated statistical quantities are. In |DF07j the 
volume preservation of the projected Heun method and of the Implicit Midpoint Rule was 
studied, by computing the determinant of the Jacobian of the discrete map, 9q"'''^/9q"'. 
It was observed that a projected Heun method (preserving the integrals), computed the 
expected average parameters, with a far too fast convergence. Further, it was observed 
that the volume was strongly contracted, implying that only a smaller region of the space 
was visited by the trajectory. They also tested thoroughly the Implicit Midpoint Rule, 
which preserved intrinsically the integrals to machine precision. The IMR is a geometric 
integrator and backward error analysis guarantees that it solves a nearby Hamiltonian 
problem for exponentially long time. The IMR does not preserve volume however it is 
symplectic, therefore it preserves the sum of 2-dimensional projected areas: for a shorter 
simulation time and a smaller number of particles (A^ = 12), and it was seen to slightly 
increase (about 3 x 10^^ for 10^ steps with stepsize 0.1). 

In our experience, if the statistical parameters were required up to d decimals, the 
relative error for the energy and enstrophy should be no larger than 10"*^ in magnitude. 
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